! a new discontinuous Galerkin method for wave equation

module WEmode
    implicit none
contains
    subroutine we_stiffmtx(pts,qdrs,m,n, p, ep, ga, at,iat,jat, a,ia, ja, am, iam, jam)
        use myfem, only: mf_localstiffmtx, mf_jacobi_qdr
        use common_functions, only: cf_lglnodes
        implicit none
        double precision, intent(in) :: pts(:,:)
        integer, intent(in) :: qdrs(:,:)
        integer, intent(in) :: m, n, p
        double precision, intent(in) :: ep, ga
        integer, intent(inout) :: iat(:), jat(:), ia(:), ja(:)
        double precision, intent(inout) :: at(:), a(:)
        integer, intent(inout), optional :: iam(:), jam(:)
        double precision, intent(inout), optional :: am(:)

        integer :: nq, dofs, ldofs, pp
        double precision :: Adump(1,1)
        double precision, allocatable :: lAe(:,:,:), lAed(:,:,:), lAedx(:,:,:), lAedy(:,:,:)
        double precision, allocatable :: xn(:), w(:)
        double precision, allocatable :: xl(:), yl(:), wl(:)
        double precision, allocatable :: G(:,:), g11(:,:), g12(:,:), g21(:,:), g22(:,:)
        integer :: ii

        nq = size(qdrs,1)
        ldofs = (p+1)**2
        dofs = ldofs * nq

        allocate(lAe(ldofs,ldofs,6), lAed(ldofs,ldofs,8))
        allocate(lAedx(ldofs,ldofs,8), lAedy(ldofs,ldofs,8))
        call mf_localstiffmtx(p,Adump,Adump,Adump,lAe,lAed,lAedx,lAedy)

        pp = p + 3
        allocate(xn(pp+1), w(pp+1))
        call cf_lglnodes(pp,xn,w)
        xn = (xn+1.d0)/2.d0
        w = w/2.d0
        allocate(xl((pp+1)**2),yl((pp+1)**2),wl((pp+1)**2))
        do ii = 1, pp+1, 1
            xl((ii-1)*(pp+1)+1:ii*(pp+1)) = xn
            yl((ii-1)*(pp+1)+1:ii*(pp+1)) = xn(ii)
            wl((ii-1)*(pp+1)+1:ii*(pp+1)) = w*w(ii)
        end do
        deallocate(xn,w)
        allocate(xn(p+1), w(p+1))
        call cf_lglnodes(p,xn,w)
        xn = (xn+1.d0)/2.d0
        w = w/2.d0

        allocate(G(nq,(pp+1)**2),g11(nq,(pp+1)**2),g12(nq,(pp+1)**2))
        allocate(g21(nq,(pp+1)**2),g22(nq,(pp+1)**2))
        call mf_jacobi_qdr(pts,qdrs,xl,yl,G,g11,g12,g21,g22)




        if ( allocated(lAe) ) deallocate(lAe)
        if ( allocated(lAed) ) deallocate(lAed)
        if ( allocated(lAedx) ) deallocate(lAedx)
        if ( allocated(lAedy) ) deallocate(lAedy)
        if ( allocated(xl) ) deallocate(xl)
        if ( allocated(yl) ) deallocate(yl)
        if ( allocated(wl) ) deallocate(wl)
        if ( allocated(G) ) deallocate(G)
        if ( allocated(g11) ) deallocate(g11)
        if ( allocated(g12) ) deallocate(g12)
        if ( allocated(g21) ) deallocate(g21)
        if ( allocated(g22) ) deallocate(g22)
    end subroutine we_stiffmtx



end module WEmode















program example_WE
    use myfem
    use common_functions
    use WEmode
    implicit none
    double precision :: O(2), width, height
    integer :: p, m, n
    double precision, allocatable :: pts(:,:)
    integer, allocatable :: qdrs(:,:)
    double precision :: ep, ga
    double precision, allocatable :: at(:), a(:), am(:)
    integer, allocatable :: iat(:), jat(:), ia(:), ja(:), iam(:), jam(:)

    integer :: ldofs, dofs

    O = (/0.d0,0.d0/)
    width = 1.d0
    height = width

    p = 1
    m = 4
    n = m
    ep = 0.d0
    ga = 1.d0


    allocate(pts((p*m+1)*(p*n+1),2), qdrs(m*n,4))
    call mf_cartesianmesh_rowelements(O,width,height,m,n, pts, qdrs)

    ldofs = (p+1)**2
    dofs = ldofs*m*n
    allocate(at(ldofs**2*m*n*2),iat(2*dofs+1), jat(ldofs**2*m*n*2))
    allocate(a(ldofs**2*m*n*2),ia(2*dofs+1), ja(ldofs**2*m*n*2))
    call we_stiffmtx(pts,qdrs,m,n, p, ep, ga, at,iat,jat, a,ia, ja)


end program example_WE
